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ABSTRACT 



Hydromagnetic stresses in accretion discs have been the subject of intense theoretical research over the past one and a half decades. 
Most of the disc simulations have assumed a small initial magnetic field and studied the turbulence that arises from the magne- 
torotational instability. However, gaseous discs in galactic nuclei and in some binary systems are likely to have significant initial 
magnetisation. Motivated by this, we performed ideal magnetohydrodynamic simulations of strongly magnetised, vertically stratified 
discs in a Keplerian potential. Our initial equilibrium configuration, which has an azimuthal magnetic field in equipartion with thermal 
pressure, is unstable to the Parker instability. This leads to the expelling of magnetic field arcs, anchored in the midplane of the disc, to 
around five scale heights from the midplane. Transition to turbulence happens primarily through magnetorotational instability in the 
resulting vertical fields, although magnetorotational shear instability in the unperturbed azimuthal field plays a significant role as well, 
especially in the midplane where buoyancy is weak. High magnetic and hydrodynamical stresses arise, yielding an effective a-value 
of around 0. 1 in our highest resolution run. Azimuthal magnetic field expelled by magnetic buoyancy from the disc is continuously 
replenished by the stretching of a radial field created as gas parcels slide in the linear gravity field along inclined magnetic field lines. 
This dynamo process, where the bending of field lines by the Parker instability leads to re-creation of the azimuthal field, implies that 
highly magnetised discs are astrophysically viable and that they have high accretion rates. 

Key words. Accretion, accretion disks - Galaxy: center - Instabilities - MHD - Turbulence 



1. Introduction 



Since the seminal work of iBalbus & Hawlevl d!991l) it has 
been widely recognised that hydromagnetic stresses play a 
central role in the dynamics of Keplerian gaseous discs. 
Much of the subsequent theoretical work has been devoted 
to the study of the magnetic dynamo driven by a combi- 
nation of magn etorotational instability (MRJ), rotation, and 
. £J | Ke plerian shear jBrandenburget aflll995t lHawlev et all Il996l 
K> see iBrandenburg & Subramania nl l2005l for an excellent review 
■ of astrophysical dynamo theory). Typically, numerical simula- 
te ' tions of the MRI-driven dynamo begin with an initial zero net 
flux magnetic field with an associated pressure which is a small 
fraction of the thermal pressure. The idea is that the MRI-driven 
turbulence should increase the characteristic coherence length 
of the magnetic field, and grow its mean strength to a signifi- 
cant fraction of the equipartition value. The simulations, how- 
ever, have had at best a mixed success in explaining high and 
persistent inflow rates observed in astrophysical accretion discs. 
Firstly, a number of recent numerical experiments have shown 
that the effectiveness of the dynamo depends in a critical way 
on the value magnetic Prandtl number, defined as the ratio of 
the collisional kinematic viscosity to the magn etic diffusivity 
dLesur & Longaret7Il2007l: l Fromang et al.Ll2007l fo l lowing sug- 
gestive non-disc simulations of Schekochih in et al and a 
careful study of pre viously published MRJ-turbulence results by 
iPessah et al.L |2007|) . The upshot of this work is that when the 
Prandtl number is significantly less than one, as is expected 



in most accretion discs dBalbus & Henril I2008I) . the dynamo 
seems to faiQ. Secondly, in the zero net flux simulations where 
all parameters are chosen so that the dynamo works, the mea- 
sured value of the Shakura-Sunyaev parameter a is typically 
of order 10~ 3 and at most 10~ 2 (e.g . IBrandenburg et all I1995 ; 
Sano et al.L [20041: iJohansen & Klahrl 120051: iFromang & Nelson . 
20061) . This is 1-2 orders of magnitude smaller than what is re- 



quired to explain t he high accretion rates in dwarf novae systems 
dKing et al.Ll2007l) . 

The weak initial magnetic fields assumed in almost all disc- 
MRI simulations may be unrealistically small for many astro- 
physical discs. We give here two examples: 
1. Extended accretion discs around central supermassive black 
holes (such as the ones traced by megamasers in nea rby galax- 
ies; see e.g. lGreenhillil20Q7Hviemmings et aUl2007l) are fed by 
molecular material in the interstellar medium. Molecular clouds 
are known to have large scale superthermal magnetic fields, 
and thus the initial magnetic fields in AGN discs are likely to 
be comparable to or larger than the thermal equipartition val- 
ues. The ~2 pc molecular circumnuclear disc in our galactic 
centre is permeated by large scale equipartition magnetic fields 
dWardle & Konigllll990trHildebrand etaTlll990t) which are cer- 
tain to play an important dynamical role in its subsequent evolu- 
tion. 



Send offprint requests to: 
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1 Also the magnitude of the turbulent magnetic stresses decrease 
as the grid resolution increases (decreasing t he numerical diffusivity 
accor dingly), with no convergence in sight dFromang & Papaloizoul . 
I2007h . In may well be, however, that this problem is related to an insuf- 
ficient size of the simulation domain, and that bigger, stratified shearing 
boxes will show better convergence dRegev & Umurhanl.[2008T) . 
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2. Accretion discs in close binary systems, which are fed by a 
Roche-lobe overflow from a tidally distorted low mass star, may 
be initially magnetised. The magnetisation 1//? of the gas as it 
departs from the donor star can be estimated as 



1 



B Snpc* 



(1) 



where B* (measured in Gauss) is the stellar magnetic field at the 
surface, c, is the speed of sound at the stellar surface, and p is the 
density of gas as it becomes detached from the star. The density 
p can be estimated as 



M 



(2) 



where M is the average accretion rate of the companion, and 
is the radius of the star. We get 



1 / M 
- ~ 0.3 5 



3km/s 



f^ff-^f (3) 
\1G/ \10 11 cm/ K ' 



The streaming and shearing motion after the gas detaches from 
the star can further amplify the magnetic field. Clearly, there is a 
realistic parameter range where the initial magnetisation is high. 

Motivated by these astrophysical considerations, in this pa- 
per we perform numerical experiments on gas discs which con- 
tain initially strong magnetic fields, with magnetic pressure com- 
parable to that of the gas. Specifically, we initialise the az imuthal 
field so that it is subject to the Parker instability (PI, iParkerl 
119661) in the vertical stratification. 

In te rms of initial conditions our p h ysical set up is 
close to IMachida. Havashi. & Matsumotol (l2000l) . However, 
iMachida et al.l d2000t) focused on the formation and evolution of 
a disc corona, and the spirit or their numerical experiment was 
different from ours. They have simulated the whole circular disc, 
and have introduced reflecting boundary conditions at the mid- 
plane of the disc, which suppresses field anchoring in the mid- 
plane. By contrast, our purpose is to understand the long-term 
dynamics of the disc, possible field confinement, and the dy- 
namo processes related to strong magnetic fields. Therefore, our 
shearing-box simulations focus on a small part of the disc and 
study it with high numerical resolution. We simulate the fluid 
both below and above the midplane, and thus have no artificial 
boundary conditions at the midplane of the disc. 

We find that both the short term and, more importantly, the 
long term behaviour of initially strongly magnetised discs is 
radically different from that of their weakly magnetised coun- 
terparts. We observe the following three-step dynamics: (a) the 
Parker instability expels azimuthal field in huge arcs, creating 
vertical field which becomes the seed for a strong magnetoro- 
tational instability, (b) matter sliding down inclined field lines 
stretches the azimuthal magnetic field and creates a vertically 
dependent large scale mean radial field, and (c) the Keplerian 
shear recreates azimuthal field from the stretching of the ra- 
dial field. The latter step close s the dynamo cycle, in much the 
same way as was sketched by iTout & Pringld d 19921) a decade 
and a half ago, although non-axisymmetric magnetorotational 
instability in the azimuthal field also plays a n important role in 
creating accretion stresses in our simulations dBalbus & Hawle 
ll992tlFogiizzo& Tag gejl 119951 : iTerquem & Papaloizoul \199 ~ 
The azimuthal field remains strongly concentrated towards the 
disc midplane; this is in contrast with the simulations of Parker 
instability whic h do not include strong Keplerian shear in the 
radial direction dKim et al.L 1 19981) . We show that the dynamo is 



robust and stable over at least tens of orbital periods, and that 
the accretion torque increases if we use a finer grid. We observe 
a- values as high as 0. 1 in our highest resolution run. 

The plan of the paper is as follows. In §|2] we detail the me- 
chanics of our numerical experiments, and in the following sec- 
tion (§|3]i we test our computer code by comparing the results 
of Parker instability simulations without Keplerian shear to the 
extensive literature that exists on the PI under rigid rotation. In 
§|4] we perform simulations with the Keplerian shear included, 
and present our main results. We devote the next section (§0 
to analysing the confinement of azimuthal flux to the disc by a 
dynamo process. In §|6]we conclude with the discussion of the 
astrophysical implications and possible future improvements of 
our work. 



2. Description of the numerical experiment 

In this section we describe the numerical method we use to solve 
the equations of ideal magnetohydrodynamics, the set up of an 
initial condition that is unstable to the Parker instability, and 
technical issues such as dissipation terms and boundary condi- 
tions. We use the Pencil CodeQ to solve the relevant partial dif- 
ferential equations. 



2.1. Dynamical equations 

We consider a local corotating patch of a Keplerian accretion 
disc. The coordinate axes are oriented such that x points radially 
away from the central gravity source, y points along the main 
orbital flow, while z points vertically out of the disc parallel to 
the Keplerian rotation vector Q. In the shearing sheet approxi- 
mation the equation of motion, for the velocity field u relative to 
the Keplerian flow, is 

— + (k ■ V)m + uf ] ^- = 2Qu y e x - (2 - q)Qu x e y - Q 2 ze z 

+—JxB--VP + f v (u). (4) 
P P 

Here q - d In Q/d In r is the shear parameter, with q = for 
rigid rotation and q = 3/2 for Keplerian rotation. The left hand 
side includes advection both by the velocity field u itself and by 
the linearised Keplerian flow mJ 0) = -qQx. The first two terms 
on the right hand side represent the Coriolis force in the x- and 
y-directions, modified in the y-component by the radial advec- 
tion of the Keplerian flow, u y = -u x duf } /dx. The third term is 
the linearised vertical component of the gravity of the central ob- 
ject, while the Lorentz and pressure gradient forces appear in the 
usual way. The high order numerical scheme of the Pencil Code 
has very little numerical dissipatio n from time-stepping the ad- 
vection term dBrandenburgj I2003I) . so we add explicit viscosity 
through the term f v (u), described in detail in j ^2.2| 

2.1.1. Alfven limiter 

In this paper we consider stratified models with 4-16 orders 
of magnitude of range in densities. In low-density regions the 
Alfven speed 



B 



VA 



(5) 



See 



http : //www . nordita . dk/data/brandenb/pencil- code/ 
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may become very high when strong magnetic fields are trans- 
ported away from regions near the midplane to form a low den- 
sity corona, and we must reduce va artificially in order to avoid 
very low time-steps in the explicit numerical scheme of the 
Pencil We have done so by replacing the density p in the 

Lorentz force term of Eq. by a modified density p* defined as 



1 _ 1 

P* P 



1 



2 \" 



M 2 



-l/n 



(6) 



Here M is a limiting value for the Alfven speed and n is an index 
that smooths the transition from the regime where the Alfven 
speed is not affected (va < M) to the region where the limiter 
applies (va > M). This yields an effective Alfven speed ~ va 
in the limit va <k M, while in the limit of high Alfven speeds the 
effective speed is 



B 



B M 



Va = 



= M for v A » M . 



(7) 



We used a smoothing index of n = 5. We are grateful to T. 
Heinemann for imp lementing this form of the Alfven limiter 
in the Pencil Code dHeinemann et all 20071). A s imilar Alfven 
speed limiter was used by iMiller & Stond (120001) . The Alfven 
limiter in the form presented here does not conserve momen- 
tum, but we have experimented with M 2 = 10 and M 2 = 100 
and found no qualitative differences in our results. Thus we ap- 
ply M 2 = 100 to all 2-D simulations, and M 2 = 10 to the 3-D 
simulations to get a longer time-step in those computationally 
expensive runs. 

2.1.2. Induction equation 

The magnetic vector potential A is evolved through the induction 
equation 



dA dA 

— + u\ 0) — = u x B + qQA y e x + f (A) . 
ot y ay 1 



(8) 



The explicit resistivity term f n (A) is discussed in $2.21 Working 
with the vector potential A rather than the magnetic field B has 
the advantage that B = V x A is kept solenoidal (V • B = 0) 
by mathematical construction. The stretching of the radial com- 
ponent of the magnetic field by the Keplerian shear is repre- 
sented by the second term on the right hand side of Eq. ([8]l 
[Brandenburg et al., 1995]. The current density /, needed for 
the Lorentz force in Eq. @, is calculated from Ampere's law 
J = jUq 1 V x B. We set the vacuum permeability po = 1. 

2.1.3. Continuity equation and equation of state 

The mass density is evolved in its logarithmic form 



<91np 
dt 



+ u ■ V lnp + u" 



(0) 



<91np 
~8y~ 



— — V ■ u + / D (lnp) . 



(9) 



The last term on the right hand side is an explicit diffusion 
term (see {32. 21 ). It is advantageous to evolve the logarithmic 
mass density when a large dynamical range in density is con- 
sidered. Also, in the initial magnetohydrostatic equilibrium, In p 
varies parabolically with z, and the Pencil Code's finite differ- 
ence approach to spatial derivatives yi elds a perfect derivative of 
parabolic curves dBrandenburgH2003h . 



3 In these magnetised, low-density regions the field is almost force- 
free, and no dynamical Alfven waves are excited. 



We close the dynamical equation system by applying an 
isothermal equation of state with P = c 2 p. The sound speed c s is 
assumed constant. 



2.2. Dissipation terms 



as 



(10) 



The viscosity term f v of Eq. appears in full generality 

f v = y 3 [V 6 H + (5 (3) ■ Vlnp)] 

+v sho ck [VV • u + (V ■ «)(Vlnp)] + Vv sho ck(V ■ «) . 

The viscosity includes both hyperviscosity with coefficient V3 
and shock viscosity with coefficient v s h oc k- A simplified third or- 
der rate-of-strain tensor S (3) is here defined as 

t-ij- <»> 

The shock viscosity coefficient is obtained by taking the negative 
part of the divergence of u, then taking the maximum over three 
grid cells in each direction, and finally smoothing over three grid 
cells, to obtain 

v s h ck = c sh0 ck<max[-V • u]+)(5xf . (12) 

Here c s h oc k is a parameter of order unity dHaugen et al.L 12004b . 

For the resistivity term / in Eq. (0 we include both hyper- 
resistivity and shock resistivity, 



773 V 6 A + 77 shock V 2 A + V77 shock V • A 



(13) 



This form is the resistivity conserves all components of the mag- 
netic field B. We set 77 shc ,ck = v s h oc k- 

Mass diffusion consists of shock diffusion, 

/ D (lnp) = D shock [V 2 lnp + (Vlnp) 2 ] + VD shock ■ Vlnp, (14) 

where £> s hock = v s h oc k is defined in Eq. ( TTZl i. We also upwind 
the finite differencing of the advection term in the continuity 
equation dDobler et al.L l2006). Mass diffusion and upwinding are 
necessary to damp out spurious small scale modes that are left 
behind due to the dispersion error in the finite differencing of the 
advection term in the continuity equation. 

We set shock diffusivity coefficients to c s h oc k = Cshock = 
^shock = 6.0 in all simulations to dissipate energy in shocks 
forming far away from the midplane of the disc. In run 
S3D_256_Lzl8, which is extended vertically compared to the 
other 3-D simulations, we had to double c s h oc k to dissipate 
enough energy in the regions above six scale heights from the 
midplane. 

We test the dependence of our results on dissipation parame- 
ters by keeping the mesh Reynolds number approximately con- 
stant with increasing resolution, i.e. scaling hyperdissipation pa- 
rameters by {6x) s and shock dissipation parameters by (Sx) 2 . 
Higher resolution simulations thus probe the effect of decreas- 
ing the dissipation coefficients. 

2.3. Initial conditions 

We assume that the initial magnetic field is toroidal, a reasonable 
assumption since Keplerian shear efficiently generates coherent 
toroidal field from a small radial component. Furthermore, we 
tune the initial field so that (a) magnetic pressure is the constant 
fraction /T 1 of the gas pressure, and (b) vertical hydrostatic equi- 
librium is enforced, i.e. 



2 \ 



= -p( z )Q 2 z --\p+-± 



dP 

= -p(z)Q 2 z-(l+p- 1 )—. (15) 
oz 
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Table 1. Simulation parameters. 



Run 


L x x Ly x L- 


N x x N y x N z 


P 


<5.:> 


Q 


1 


Af 


R2D.256 


0.0 x 24.0 x 12.0 


1 X 128 x 256 


1.0 


0.0 


1.0 





100 


R2D.512 


0.0 x 24.0 x 12.0 


1 x 256 x 512 


1.0 


0.0 


1.0 





100 


R2D_256Xz24 


0.0 x 24.0 x 24.0 


1 x 256 x512 


1.0 


0.0 


1.0 





100 


R3D.256 


12.0 x 24.0 x 12.0 


128 x 128 x 256 


1.0 


0.0 


1.0 





30 


S3D.256 


12.0 x 24.0 x 12.0 


128 x 256 x 128 


1.0 


0.0 


1.0 


3/2 


30 


S3D.512 


12.0 x 24.0 x 12.0 


256x 512x256 


1.0 


0.0 


1.0 


3/2 


20 


S3D_256_b3 


12.0 x 24.0 x 12.0 


128 x 256 x 128 


3.0 


0.0 


1.0 


3/2 


30 


S3D_256_Lzl8 


12.0 x 24.0 x 18.0 


128 x 256 x 192 


1.0 


0.0 


1.0 


3/2 


20 


S3D_256_Bz0.03 


12.0 x 24.0 x 18.0 


128 x 256 x 192 


1.0 


0.03 


1.0 


3/2 


20 



The first column gives the name of the simulation, while the box size, in units of the thermal scale height H = c s Q~' , and the grid resolution 
are given in the following two columns. The parameters given in the last five columns are: the ratio of thermal to magnetic pressure /? in the 
initial azimuthal field, the mean imposed vertical field (B z ) in units of the midplane azimuthal field, the angular frequency Q of the box, the shear 
parameter q, and finally the simulation time in orbits. 



Mathematically, the initial density, magnetic field, and vector po- 
tential are given by 

p(z) = p exp[-z 2 /(2//„ 2 )], (16) 



B 

= /T'c^oexpt-z 2 /^ 2 )], 
2^o p 



(17) 
(18) 



A x (z) = ^2nn Q B-^cl P(s Hpzrf[zl{2Hp)} . 

Here Hp — y/l + B^c^/Q is the gas scale height, and po is the 
initial density in the midplane. The thermal pressure alone would 
give rise to a scale height of H = H^, = c s /Q. We will use this 
more familiar scale height as our unit of length, setting H = 1 
throughout the paper. 

Since we shall consider many (six to twelve) scale heights 
above and below the midplane, we run into the problem that 
the finite differencing of the magnetic pressure gradient under- 
estimates its actual value and leads to spurious accelerations far 
away from the midplane. This is not a problem for the equilib- 
rium density stratification, since the logarithmic density varies 
as a parabola with height over the midplane, and any parabolic 
shape has a perfect numerical derivative in the symmetric finite 
difference scheme of the Pencil Code. In order to give mag- 
netic stratification equal opportunity as the density, and to ob- 
tain a perfect initial magnetohydrostatic equilibrium, we mea- 
sure magnetic field and current density relative to their initial 
values, as explained below. 

We split the field and the current into two components, as 
follows. For the Lorentz force in Eq. © and the induction term 
in Eq. © we replace 

B -> B (0) + B\P(z)e y , (19) 



J {0) + J?\z)e x . 



(20) 



The initial magnetic field B x (z) is given by Eq. ( 1171 . while the 
initial current density jf \z) is calculated as 



J?\z) = — ^ 



1 8By 

HQ dz 



splitting does not introduce any extra spurious effects. We have 
checked that the Parker instability develops similarly in mod- 
els with and without this field splitting, and have found that the 
models without splitting indeed converge towards the models 
with splitting as the resolution increases. 

We perturb the initial Keplerian velocity field by Gaussian 
noise of amplitude Su ~ 10~ 3 . This perturbation seeds the lin- 
ear Parker instability, and the noise also contains leading waves 
that can be amplified by non-axisymmetric magnetorotational 
instability in the azimuthal magnetic field in simulations with 
Keplerian shear. 

2.4. Boundary conditions 

We set the usual shearing sheet boundary conditions in the x 
and y directions (shear-periodic in x and periodic in y). At the 
upper and lower boundaries we impose free-slip conditions, 
du x /dz - du v /dz = u z = 0, for the velocity field u and per- 
fect conductor, A x = A y = dA,/dz = 0, for the magnetic vec- 
tor potential A. This effectively allows the radial and azimuthal 
components of the magnetic field to evolve freely at the upper 
and lower boundary planes, while the vertical component of the 
magnetic field is forced to be zero. Globally the average value of 
all components of the magnetic field are then conserved. Perfect 
conductor boundary conditions do not allow azimuthal flux to 
flow out of the box, and therefore there is a concern that the 
presence of the boundaries will artificially enhance the field in 
the disc. To address this concern, we have made sure that our re- 
sults have converged with increasing vertical extent of the box; 
see ^5]and Figs. 171 and [121 

For the mass density In p we set a symmetric boundary con- 
dition with din p/dz = at the upper and lower boundaries. 
Although this precludes the gas pressure gradient from partic- 
ipating in magnetohydrostatic equilibrium there, the free-slip 
boundary condition for the velocity field, with zero vertical ve- 
locity component, means that this does not cause any unwanted 
accelerations. 



^BMpdzlilHl)] exp[-z 2 /(4^ 2 )] . (21) 2 5 ' Simulation parameters 



The effect of this splitting is that the initial, mean field is not 
exposed to any (numerical or explicit) resistivity. Since we wish 
for resistive effects to go to zero with increasing resolutiorfl, the 

4 This is because our numerical resistivity is always orders of magni- 
tude greater than in realistic accretion discs. 



We perform a series of simulations of the evolution of Parker and 
magnetorotational instabilities, usually starting from an equipar- 
tition (J3 = 1) accretion disc in magnetohydrostatic equilib- 
rium (as described in jj2.3l l. Simulation parameters are written 
in Table Q] The main runs are the shearing sheet simulations 
S3D_*, representing varying resolution (S3D_256, S3DJ512), 
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Fig. 1. Evolution of the Parker instability in 2-D rigid rotation (run R2D_256). Overlaid on the density are magnetic field streamlines 
(white lines) and velocity field vectors (white arrows, averaged over 8 grid points, the longest arrows represent approximately four 
times the sound speed). The initial stratification is unstable to magnetic buoyancy, and magnetic field arcs begin to rise from the 
midplane. The arcs merge to form longer arcs, and eventually the system settles down into a new equilibrium state with two superarcs 
and four dense pockets of matter in the midplane. 



higher initial ratio of gas to magnetic pressure (S3D_256_b3), 
increased vertical extent of the box (S3D_256Xzl8) and an ex- 
ploration of the effect of a moderate amount of net vertical field 
(S3DJ256_Bz0.03). For testing purposes we also run a series 
of 2-D and 3-D rigid rotation simulations (R2D_256, R2D_512, 
R2D_256_Lz24, R3D_256), which we present in the next section. 



3. Testing the code: evolution of the Parker 
instability under rigid rotation 

In this section we study the evolution of the "pure " Parker insta- 
bility in 2-D and 3-D simulations without shear dParkerlTl966t 
iKim et aflll997b . This allows us to test our code against the pre- 
viously published results, and to review previous work done on 
the evolution of the Parker instability in systems with no system- 
atic shear. 

iBasu et aT] d 1997b performed 2-D inertial frame simulations 
of the non-linear evolution of the Parker instability. They found 
that primarily antisymmetric midplane crossing modes are ex- 
cited by the random noise initial condition. The instability devel- 
ops quickest in the tenuous gas high over the midplane. As the 
magnetic field lines rise in big arcs, matter streams unobstructed 
down the lines that are no longer horizontal. Dense pockets of 
gas (approximately a factor two times more dense than the av- 
erage state) form in the midplane, anchoring the rising mag- 
netic field lines there. The final state of the Parker instability, 
in which the tension in the bent magnetic field lines balances the 
buoyancy, has a low er energy than the initial equilibrium state 
dMouschoviasl 1 1 9741) and is therefore preferred. 

In Fig. [T] we show the evolution of the Parker instability ob- 
tained with the Pencil Code under rigid rotation in the 2-D (y, z) 
plane (run R2DJ256). Magnetic field lines initially rise in arcs 



3.0 L 

2.5 : \ 

i.o ■ — - 

0.5 '- 256x128 : 

: 512x256 : 

nn f ■ 

20 40 60 80 100 

t 

Fig. 2. The maximum gas density as a function of time for 2- 
D simulations of the Parker instability (under rigid rotation). 
Increasing resolution leads to no significant change of the max- 
imum density, showing that the Parker instability operates in a 
robust way already at 256 x 128 grid points. 

of wavelength A-pi as 1H. As matter streams sideways along the 
rising field lines (see arrows in Fig. [TJ dense clumps of gas ap- 
pear in the midplane after 10 orbits. Eventually the magnetic arcs 
merge and form two superarcs with four dense gas clumps in the 
midplane. 
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2-D rigid rotation 3-D rigid rotation 




Fig. 3. Mean azimuthal magnetic field, averaged over the y- 
direction, as a function of height z over the midplane in 2-D 
simulations with rigid rotation. The initial Gaussian magnetic 
pressure profile is lifted by the Parker instability, reducing the 
azimuthal flux by a factor of approximately 1 .5 in the midplane. 
The flux is transported to above four scale heights from the mid- 
plane, and an approximately force free corona develops at these 
high altitudes. Increasing the vertical extent of the computational 
domain from 12 to 24 scale heights causes the flux to be trans- 
ported slightly further away from the midplane, but the overall 
shape of the magnetic field distribution stays approximately the 
same. 



Fig. 4. Mean azimuthal magnetic field versus height over the 
midplane for the Parker instability in 3-D rigid rotation. The 
azimuthal field quickly spreads evenly over the entire vertical 
extent of the box, in contrast to the 2-D case (shown in Fig. [3]) 
which makes a non-linear transition to a state that still has az- 
imuthal flux differences. 



azimuthal magnetic field distribution is nevertheless relatively 
unchanged. 



In Fig. [2] we show the maximum gas density in the box as 
a function of time. The clumps in the midplane are overdense 
by around factor 1.6 - 2.2 compared to the initial state, and this 
contrast is relatively well converged already at 256 x 128 grid 
points (as is the non-linear evolution of the Parker instability in 
ge neral). Our simulat ions are in good qualitative agreement with 
the lBasu et al .1(1 19971) simulations. 

3.1. Evolution of the azimuthal field 

In Fig.[3]we show the mean azimuthal magnetic field as a func- 
tion of height over the midplane in the saturated state of the 2- 
D Parker instability under rigid rotation. The initial azimuthal 
field is redistributed by the Parker instability so that some of the 
midplane flux ends above four scale heights from the midplane. 
This scale is similar to the typical azimuthal wavelength of the 
Parker instability, so the final scale of vertical variation of the 
field likely comes about as the magnetic tension in the rising 
field lines grows to match the buoyant rise of the field lines. The 
non-linear transition from the initially stratified state to a state 
with magnetic arcs anchored in the midplane leads to the forma- 
tion of an approximately force-free corona extending from ~4 
scale heights from the midplane. In the corona, the magnetic ten- 
sion and pressure terms are large (compared to their difference) 
and oppositely directed, and current flows parallel to the mag- 
netic field lines. Increasing the vertical extent of the box from 
12 to 24 scale heights leads to a transport of some magnetic flux 
to regions with \z\ > 6H, at the cost of a 20-30% decrease in the 
magnetic flux in the region |z| < 6H. The overall shape of the 



3.2. Rigid rotation in 3-D 

Two-dimensional models of the Parker instability allow only the 
undular mode to develop , but not the related interchange mode 
dMatsumoto et al.L 1 1993b . and hence the field has no possibil- 
ity to escape from the disc. In three-dimensional rigidly rotating 
disc simulations, the field does escape. In Fig. [4] we show the 
evolution of the azimuthal field in a 3-D model including rota- 
tion (but still no shear). Now the azimuthal field quickly redis- 
tributes itself evenly with height over the midplane, leading to a 
state where vertical gravity is balanced only by thermal pressure. 
Allowing the mixed undular/interchan ge instabili t y to ta ke place 
leads to a genuinely turbulent state dKim et all 1 1998b . which 
mixes the azimuthal field vertically and eventually spreads it 
uniformly throughout the entire extent of the box. In the next 
section we will show, however, that Keplerian shear completely 
suppresses the vertical spreading of the azimuthal field, which 
remains concentrated towards the disc midplane due to a dynamo 
coupling between B x and B y . 

The Parker instability in a shearing, rotating frame (which is 
the relevant frame for most astrophysical purposes) must neces- 
sarily be more complicated than in the inertial frame. The intrin- 
sically non-axisymmetric instability is deformed by the radial 
shear, while the whole coordinate frame can tap energy from the 
gravitati onal potent ial through Maxwell and Reynolds stresses. 
Already IShul d 19741) expanded the Parker instability to include 
the effect of shear. The conclusion of the analytical stability anal- 
ysis is that no amount of shear can completely stabilise against 
the Parker instability, as the instability can develop at such small 
radial scales that growth happens faster than shear. 
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Fig. 5. The vertical magnetic field at the sides of the simulation box for run S3D_512. The initial fastest growing wavelength 
of the Parker instability is around 5 scale heights in the azimuthal direction and 0.5 scale heights in the radial direction. The 
magnetorotational instability eventually sets in as the vertical field grows in amplitude, and the box turns turbulent with a strongly 
fluctuating vertical field. 



4. The Parker instability in Keplerian shear 

Having gained insight into the evolution of the Parker instability 
in 2-D and 3-D models with rotation but no shear, we now turn to 
the combined effect of Parker and magnetorotational instabilities 
in a shearing sheet model of a Keplerian accretion disc. 



4.1. Creation of vertical field 

In Fig. [5] we show the evolution of the vertical component of 
the magnetic field at the sides of the box for run S3D_512. 
The following sequence of events is observed: first the (sheared) 
Parker instability grows unaffected by the MRI, but as the verti- 
cal magnetic fields gets a significant amplitude, the MRI sets in 
and leads to a turbulent state where both PI and MRI determine 
the dynamics. The Parker instability initially sets in high up in 
the atmosphere where the buoyancy is strongest, at the typical 
wavelength A y = 5H and X x = 0.5//. We actually observe that 
the radial wavelength of the Parker instability decreases with 
increasing resolution as the reduced dissipation allows growth 
at higher wave numbers. Eventually the whole disc goes into a 
non-linear turbulent state with strongly fluctuating vertical fields 
and twisted magnetic field loops extending far above the mid- 
plane of the disc. The typical radial coherence scale in this state 
is much longer than in the initial linear growth phase (compare 
t = 2.0T olb to t = 15.0T OTb in Fig.|5}. 



The initial azimuthal field 
sient) non-axisymmet ric 
itv" dBalbus & Hawlevi 
iTerquem & Papaloizoul 



is also unstable 



1992; 



1996 



ma gnetorotational 
iFpglizzo & Taggerl 



to a (tran- 
"instabil- 



Keppens et all 



growth rate for modes with a small vertical wavelength^- The 
absence of high k z modes in the first panel of Fig. [5] indicates 
that non-axisymmetric MRI is not the dominant driver of the 
linear dynamics of the system. However, in Fig.|6]we show the 
vertical velocity field at the y = -L v /2 plane for run S3D_512 
during the linear growth phase. The top panel exhibits the usual 
signs of Parker instability high up in the atmosphere, but chang- 
ing the contrast by two orders of magnitude (bottom panel) 
reveals modes of very short radial and vertical wavelength in 
the midplane of the disc, representing trailing waves that have 
already been amplified by magnetorotational instability in the 
azimuthal field. Transient amplification of swinging waves may 
also be important in deter mining the non-line a r evolution of the 
system at later times (e.g. lRincon et all l2007t iLesur & QgilvieL 
120081) . 

Although production of vertical field is not a main ingredient 
of azimuthal MRI, any unstable mode with a non-axisymmetric 
vertical velocity component can create significant vertical field 
by stretching the large scale azimuthal field (although highly 
variable in the vertical direction). The production of large scale 
vertical field by the Parker instability, on the other hand, is 
demonstrated in Figure [7] Here we plot the mean azimuthal 
field and the root-mean-square of the vertical field as a func- 
tion of height over the midplane for the mixed PI/MRI shearing 
sheet simulations. Significant vertical magnetic fields of magni- 



Magnetorotational instability in the azimuthal field has highest 



5 It was shown by iFoglizzo & Tagged d!994 119951) that in the p = 
1 case the azimuthal field MRI should have a slightly higher growth 
rate than the Parker instability. However, their linear stability analysis 
assumed a constant gravity field, rather than a linea r field which is the 
more relevant for Keplerian discs dKim et all 1 19971) . Hence the growth 
rate of the Parker instabi lity several scale heights from the midplane 
may be underestimated in lFoglizzo & Tag geH dl994l) . 
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Fig. 7. Confinement of the azimuthal field in shearing box sim- 
ulations. The mean azimuthal field and the root-mean-square of 
the vertical magnetic field are plotted as a function of height 
over the midplane. The Parker instability creates strong vertical 
fields from the initial azimuthal flux (dotted line in top plot). 
Vertical loss of azimuthal flux is nevertheless completely sup- 
pressed outside of a field reversal occurring at approximately 
four scale heights from the midplane. Doubling the resolution 
(S3D_512) leads to an increased vertical field strength and a de- 
creased azimuthal field in the midplane, due to a significant in- 
crease in turbulent motions at higher resolution (see Table 0. 



Fig. 6. The vertical velocity component at the y = -L y /2 plane 
during the linear phase of run S3D_512 (at t = 2.4r ol -b). The 
top panel shows the high radial wavenumber structures typical 
for the 3-D Parker instability. However, changing the contrast 
by two orders of magnitude reveals structures of high vertical 
wavenumber in the midplane of the disc, likely a result of mag- 
netorotational instability in the azimuthal field. 



tude fi-~0.2 (with associated pressure B Z ~5Q) arise due to mag- 
netic buoyancy. There is a reversal of the azimuthal magnetic 
field at around four scale heights from the midplane. The cor- 
respondence between run S3D_256 and run S3D_256Xzl8 is 
extremely close, so the sign reversal is not due to the presence 
of the vertical boundary. Increasing resolution leads to some in- 
crease in vertical field strength in the midplane, but the field 
reversal still takes place at the same height over the midplane. 
Although the non-linear state of the combined PI and MRI is 
time-dependent and fluctuating, there is a clear saturation and 
vertical confinement of the mean azimuthal magnetic field (we 
return to the issue of field confinement in $5[. 



4.2. Stresses 

The magnetorotational instability in turn feeds off both the verti- 
cal field lines created by the Parker instability and the azimuthal 
field component still present in the non-linear state. The result- 
ing Maxwell stress (B X B V ), averaged over time and radial and 
azimuthal directions, is shown as a function of height over the 



midplane in Figure [8] We have divided the stress into contribu- 
tions from the mean field, ((B x ) X y(By) xy )i, and fluctuating field, 
((B x - (B x )xy)(By - (B v ) xv )) xyt . The occurrence of a large scale 
B x is discussed in detail in §5\ The radial field couples with the 
large scale azimuthal field to yield a large scale structure in the 
Maxwell stress as well, relatively well converged when increas- 
ing the resolution. 

The fluctuating magnetic field is associated with significant 
stresses around the midplane, with a Maxwell stress in the range 
0.03 . . .0.1 in the regions within a couple of scale heights from 
the midplane. The stress from this fluctuating field increases sig- 
nificantly with increasing resolution, as the decreasing dissipa- 
tion on small radial scales allows the Parker instability to create 
stronger vertical fields, while the decreased numerical and arti- 
ficial dissipation lets the both vertical field and azimuthal field 
MRI develop faster. 

In Fig.|9]we show the measured mean Maxwell stress (B X B Y ) 
as a function of time. The run S256_3D_Bz0.03 has a moder- 
ate vertical field (B z = 0.03) imposed through the bo:xQ. With 
this set up the magnetorotational instability sets in before the 
Parker instability does, creating significant stresses already af- 
ter a few orbits. Eventually however the Parker instability de- 
velops as usual, and the stresses reach saturation at a level that 
is only slightly higher (in absolute value) than for the zero net 
vertical flux run. Thus it seems that it is not very important 
whether the magnetorotational instability develops before the 



6 We exclude the vertical field from the top and bottom six grid zones, 
to avoid conflict with the zero-vertical-field boundary condition. 
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Fig. 8. The Maxwell stress {B X B X ) as a function of height over 
the midplane. The top panel shows the stress from the mean 
field {B x ){By), while the bottom panel shows the contribution 
from the fluctuating field B' — B — (B). The magnetorotational 
instability feeds off both the azimuthal fields and the large scale 
vertical fields created by the Parker instability, to create signif- 
icant stresses, up to {B X B V ) * -0.1 in the midplane of the high 
resolution run S3D_512. 



Parker instability or vice versa. This situation may nevertheless 
change when going to either weaker azimuthal fields or stronger 
vertical fields, in which case the turbulent diffusion created by 
the magnetorotational instability may reduce the midplane az- 
imuthal flux quick er than the Parker instability can grow (but see 
iBlaes et ail, 120071 where the Parker instability arises in simula- 
tions with a net azimuthal field that is much weaker than ours). 
In the case of a strong, imposed vertical field, there is however 
already an in inexhaustible source of accretion stresses present 
in the disc without the need to invoke an additional mechanism 
based on azimuthal fields and Parker instability. 

In Table|2]we show the measured statistical properties of the 
non-linear state of the combined Parker and magnetorotational 
instabilities. We divide into the midplane regions with |z| < 2 
(containing 85% of the gas mass) and atmosphere regions with 
|z| > 2 (containing 15% of the gas mass). The midplane regions 
are characterised by strong magnetic fields and high accretion 
through magnetic stresses, whereas the atmosphere has stronger 
velocity fields and higher Reynolds stress, but weaker mag- 
netic fields and Maxwell stresses. Decreasing the initial mag- 
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Fig. 9. The mean Maxwell stress as a function of time. The high 
resolution run S3D_5 12 has faster growth of the stress and higher 
saturation level. Imposing the gas to a moderate net vertical field 
(dash-dotted line) lets the vertical field MRI develop from the 
beginning, but eventually the stresses saturate at an only slightly 
higher level (in its absolute value) than in the zero net vertical 
flux case. 



netic pressure by a factor of three (run S3D_256_b3) leads to an 
expected decrease in both magnetic energies and stresses. 

The z-dependent Reynolds and Maxwell stresses can 
be translated into an a verage turbulent viscosity (following 
iBrandenburg et all 1 1995b . 



ipU X Uy) 



(kin) 



<P>, 



1 

/'() 



-{B x B y ) = "-Qv ( r g \p)- 



(22) 



(23) 



Here v, •* and yj mag) are the turbulent viscosities due to the ve- 
locity field and the magnetic field, respectively. We can further 
normalise the t urbulent viscosities by the sound speed c s a nd gas 



normalise tne turbulent viscosities by tne sound speed c s an< 
scale height H (IShakura & Sunyaevlll973l:lPringlelll98lh . 



(kin) 
(mag) 



(kin) 



(24) 
(25) 



For the gas scale height we use the expression H — Hp — 

Vl +B~ l c s /Q, with B = P/P mag given by the initial value. The 
resulting a-values are given in Table [3] The large scale ver- 
tical fields arising from the Parker instability, together with a 
non-axisymmetric instability in the azimuthal field, give rise to 
a t as 0.05 . . . 0.1, resulting in very high accretion rates through 
M = 3nv,X^ JPringlelll98ll) . 

lMachida"et al. (2000) observed similarly high a- values in 
their global simulation of a strongly magnetised accretion disc, 
and they attributed the Maxwell stresses to a non-axi symmetric 
magne torotational instability of the azimuthal field. iKim et aD 
d2002l) considered the competition between Parker and grav- 
itational instabilities in the context of a galactic potential. 
Simulations including rotation and shear nevertheless did not 



10 Johansen and Levin: High accretion rates in magnetised Keplerian discs 

Table 2. Statistical flow properties in the saturated, turbulent state. 

|g| < 2H (85% of the mass) 



Run rms(u,) rms(w v ) rms(«,) mean^Wy) rms(£ v ) rms(B v ) rms(Bj) mean(B x B y ) 



S3D.256 


0.346 


+ 


0.020 


0.261 ± 


0.020 


0.344 


+ 


0.016 


0.017 ± 


0.005 


0.189 


+ 


0.009 


1.227 


+ 


0.028 


0.172 


+ 


0.005 


-0.055 


+ 


0.006 


S3DJ12 


0.366 


+ 


0.017 


0.331 ± 


0.009 


0.338 


+ 


0.025 


0.031 ± 


0.003 


0.287 


+ 


0.015 


1.143 


+ 


0.027 


0.230 


+ 


0.012 


-0.087 


+ 


0.016 


S3D_256.b3 


0.341 


+ 


0.017 


0.232 ± 


0.009 


0.359 


+ 


0.021 


0.011 ± 


0.001 


0.142 


+ 


0.009 


0.801 


+ 


0.016 


0.128 


+ 


0.006 


-0.040 


+ 


0.005 


S3D_256_Lzl8 


0.400 


+ 


0.028 


0.314 ± 


0.013 


0.370 


+ 


0.030 


0.025 ± 


0.003 


0.208 


+ 


0.011 


1.240 


+ 


0.043 


0.187 


+ 


0.013 


-0.058 


+ 


0.010 


S3D_256_Bz0.03 


0.395 


+ 


0.024 


0.292 ± 


0.013 


0.408 


+ 


0.012 


0.021 ± 


0.002 


0.216 


± 


0.011 


1.347 


+ 


0.018 


0.207 


+ 


0.003 


-0.084 


+ 


0.012 



\z\>2H (15% of the mass) 

Run rms(n r ) rms(« v ) rms(n-) mean(t( v M v ) rms(B c ) rms(B T ) rms(B-) mean(B. t B y ) 



S3D.256 


1.108 


± 


0.084 


0.696 


+ 


0.042 


0.635 


+ 


0.016 


0.059 


+ 


0.066 


0.116 


± 


0.007 


0.486 


+ 


0.015 


0.100 


± 


0.005 


-0.016 


+ 


0.004 


S3D.512 


1.072 


+ 


0.086 


0.793 


+ 


0.036 


0.665 


+ 


0.035 


0.096 


+ 


0.026 


0.163 


+ 


0.004 


0.547 


± 


0.016 


0.138 


+ 


0.006 


-0.025 


+ 


0.004 


S3D_256.b3 


1.159 


+ 


0.125 


0.634 


+ 


0.075 


0.597 


+ 


0.067 


0.022 


+ 


0.065 


0.059 


+ 


0.005 


0.194 


+ 


0.010 


0.044 


+ 


0.003 


-0.004 


+ 


0.001 


S3D_256Xzl8 


1.439 


+ 


0.092 


0.956 


+ 


0.029 


0.897 


+ 


0.077 


0.130 


+ 


0.054 


0.100 


+ 


0.005 


0.390 


+ 


0.017 


0.085 


+ 


0.006 


-0.012 


+ 


0.002 


S3DJ256-Bz0.03 


1.403 


+ 


0.101 


0.832 


± 


0.058 


0.735 


+ 


0.033 


0.192 


+ 


0.099 


0.135 


+ 


0.009 


0.448 


± 


0.027 


0.107 


+ 


0.006 


-0.016 


+ 


0.004 



The two tables indicate the root-mean-square (and its temporal fluctuation interval) of the velocity field and the magnetic field, and the mean 
product of the radial and azimuthal components of velocity field and magnetic field. The top table shows the values within two scale heights of the 
midplane, while the bottom table shows the values measured above two scale heights. 



show any sign of the magnetorotational instability and the au- 
thors note that the magnetorotational instability in the azimuthal 
field grows too slowly to show up in the simulation time-scale 
of 4-5 orbits. In this work we find that the onset of turbulence 
fits a two layer scenario where the atmosphere is dominated by 
Parker instability and MRI in the ensuing vertical fields, whereas 
azimuthal field MRI drives the linear growth within a couple of 
scale heights of the disc midplane. 

4.3. Energy spectra 

The kinetic and magnetic energy spectra are shown in Fig.fTOlfor 
runs S3D_256 and S3D_512. We have taken the power at scale k, 
P(k), and summed over shells of constant wave number k — \k\, 
excluding the anisotropic scale k = 2n/(24H) that is only present 
in the y-direction. Fig. [TUl shows that kinetic energy dominates 
over magnetic energy by approximately a factor of two at most 
scales, except for the 1-2 largest scales which are dominated by 
the magnetic field. 



10' 




k/(2%H~ l ) 



5. Field confinement 

An important issue related to our proposed path to accretion is 
whether the original azimuthal flux can stay confined in the disc 
or whether it will escape to infinity by the action of turbulent 
resistivity (i.e. magnetic buoyancy). The combination of Parker 
and interchange instabilities in the non-shearing frame will even- 
tually redistribute the azimuthal magnetic field evenly over the 
entire box height (see Fig . |4j) . We emphasise that we do not see 
the same behaviour in the shearing sheet. The magnetic field and 
velocity field stay statistically constant for the entire duration 
of the simulations, with no sign of decay or gradual loss of az- 
imuthal flux (Fig. [T|i. 

5.1. Vertical structure of the field 

In Fig.QT]we dissect the equilibrium state of B v by averaging dif- 
ferent terms from the induction equation over x and y and over 
evenly spaced snapshots between t = 10r or b and t = 30r oi b. 
There is an almost perfect equilibrium between turbulent trans- 
port of magnetic fields [B y = -V • (uB v ), representing compres- 
sion and advection due to both the mean part of the velocity field 
and the fluctuating part ("turbulent resistivity")] and azimuthal 



Fig. 10. Kinetic and magnetic energy spectrum for runs S3D_256 
and S3D_512. The power P(k) has been averaged over snapshots 
between 15 and 20 orbits and summed over shells of constant 
wave number k = \k\. Kinetic energy dominates over magnetic 
energy at most scales, except for the few largest scale of the box. 



field created by the shear from the radial field [B y = —(3/2)QB x ]. 
This w-dynamo allows for a closure of the entire process by 
which strong accretion occurs from initially purely azimuthal 
fields. The magnetic stretching term (B ■ V)m v (which does not 
include stretching by the Keplerian shear because we measure 
velocities relative to the main shear) is generally of the same 
sign as the turbulent transport term and thus acts oppositely of 
the shear term. 

The vertical and temporal dependence of the radial magnetic 
field B x is shown in Fig. [12] as a function of height over the mid- 
plane z and time t measured in orbits. After around 5 orbits a 
clear structure evolves out of the flow, with a peak of B x in the 
midplane of the disc, followed by a valley and an additional peak 
at each side of the midplane. This structure stays statistically un- 
changed for the entire duration of the run. A similar steady state 
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Run 




(km) 

y t 




ff (mag) 


a {km) 


of* 


S3D.256 


0.075 


0.020 


0.095 


0.053 


0.014 


0.067 


S3D_512 


0.107 


0.033 


0.141 


0.076 


0.023 


0.100 


S3D_256.b3 


0.042 


0.010 


0.052 


0.030 


0.007 


0.037 


S3D_256.Lzl8 


0.078 


0.017 


0.095 


0.055 


0.012 


0.067 


S3D_256.Bz0.03 


0.096 


0.023 


0.120 


0.068 


0.017 


0.085 



Table 3. Turbulent viscosity coefficients and cr-values, based on Eqs. d22l i-(l25ll. 
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Fig. 11. Different terms in the azimuthal component of the in- 
duction equation averaged over x and y and over evenly spaced 
snapshots between t = 10r oi b and t = 30r or b, for run S3D_256. 
The azimuthal magnetic field is in equilibrium between turbu- 
lent transport, representing compression and advection due to 
both the mean and the fluctuating velocity field (turbulent resis- 
tivity), and the Keplerian shear w-dynamo that creates B y out of 
B r . 



structure of the radial field was observed by lHanasz et al 1 (120021) 
in rigid rotation si mulations of the gala ctic Parker instability. In 
the simulations of lHanasz et al.1 (|2002) B x would be either pos- 
itive or negative in the midplane, depending on the radial wave 
number of the initial perturbation. In the shearing sheet the radial 
wave number of the PI is forced by the shear to be quite high, and 
our observation of a positive radial field in the mi dplane agrees 
with th e high radial wave number simulations of lHanasz et alj 
(120021) . Because we include shear in our simulations we addi- 
tionally see the constant creation of a z-dependent B v , by the 
Keplerian shear, which balances out the turbulent transport and 
prevents the azimuthal field from spreading evenly over the box 
height. 

In Fig. [13] we show the power contribution to the large scale 
B x as a function of time. We have first averaged B x and the ad- 
vection, compression and stretching terms of the radial com- 
ponent of the induction equation over the x- and y-directions. 
The power contribution of each term is extracted by multiply- 
ing the time derivative B x with B x and averaging over z. Power 
is provided almost exclusively by the magnetic stretching term 
B x — B ■ Vm v , whereas both advection and compression extract 




Fig. 12. The radial component of the magnetic field, B x , av- 
eraged over the radial and azimuthal directions, for runs with 
L z = 12H (top panel) and L : = 18// (bottom panel). The ver- 
tical structure quickly develops a peak in the midplane and one 
additional peak on each side of the midplane, and this structure 
stays statistically unchanged for the duration of the simulation. 
This radial field drives an w-dynamo that balances the turbulent 
resistivity acting on the mean azimuthal field (see Fig. ITTt. 



energy from B x at all times (except for a few peaks to around 
zero power). The magnetic stretching term provides a direct cou- 
pling between B x and B y , indicating that B x is created from B y 
in a dynamo process. 



5.2. Dynamo 

The creation of a systematic radial field component can be un- 
derstood from gas that plunges down field lines that have been 
bent by the Parker instability dParkerl 119921 lHanasz & LeschL 
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Fig. 14. Derivation of the dynamo a&yn from a single snapshot (top panels) and from the average of several snapshots (bottom panels). 
The (negative) integral of the magnetic stretching term - f B ■ Vu x dz is shown together with the azimuthal magnetic field B y , both 
as functions of height over the midplane, in the left panels. The electromotoric force contribution from the stretching term is of 
opposite sign on each side of the midplane. The right panels show the correlation between B y and - J B ■ Vu^dz, with orange (grey) 
denoting points above the midplane and blue (dark) denoting points below. There is a very clear correlation (anticorrelation) above 
(below) the midplane, indicating a dynamo a that is positive above the midplane and negative below, and of order a^yn ~ 0.023c s . 



119931: lHanasz et all |2004|) . As gas slides azimuthally along the 
field lines, the Coriolis force causes a counter clockwise rotation 
around dense clumps and clockwise rotation around underdense 
regions, twisting magnetic field lines in such a way that the per- 
turbed electromotoric force points parallel to the (unperturbed) 
azimuthal field line. Radial field is subsequently created from the 
stretching of the perturbed field lines, in much the same way as 
imagi ned for the canonical ff-mechanism dParkerl [l 955b iMoffattl 
119781) . 

To explore this scenario in more quantitative terms we follow 
IMoffattl (11978!) and expand the velocity and magnetic fields in 
constant and fluctuating parts, u-U + u',B-B + b', leading 
to the following evolution equation for the mean magnetic field 

B = Vx(UxB + a dyn B - T]J) . (26) 

Here overlines denote azimuthal and radial (and possible tempo- 
ral) averaging, while ad yn is a parameter that describes the pro- 
portionality between mean field B and fluctuating electromotoric 
force u' x V . The parameter represents turbulent resistivity, 



proportional to the current density / = V x B of the mean field. 
Writing out the radial component of Eq. d26*l i we get 

d d - d 2 B y 

B x = -y(U z B x ) - -(a iyn B y ) + nt-jl . (27) 

The first term is due to advection of the large scale radial field by 
the large scale vertical velocity component. Fig. Qj] showed that 
the source of magnetic energy in the x-component is the mag- 
netic stretching term. We may thus identify the positive contri- 
bution to B x with the properly averaged stretching term, 



B x = B Vu x . (28) 
We can subsequently estimate ttdyn from 



a dyn B y = - J B ■ V«,dz (29) 

by integrating the magnetic stretching term over z. This is of 
course only a crude approximation of that ignores many of 
the complications of analysing the evolution of the mean field 
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Fig. 13. The power contribution to the radially and azimuthally 
averaged radial field, as a function of time measured in orbits. 
Power is provided almost exclusively by the magnetic stretching 
term B ■ Vm x , while advection and compression are both sink 
terms. 



comp onent (see e.g. iBrandenburgj 1200 it iBrandenburg et all 
120081) . but this method gives a good order of magnitude esti- 
mate of the efficiency of creating radial field from the large scale 
azimuthal field. 

In the left panels of Fig.[l4]we plot the integral - J B ■ Vu x dz 
as a function of height over the midplane and compare it to 
the azimuthal magnetic field B Y . The electromotoric force con- 
tribution of the magnetic stretching term is of opposite sign 
at each side of the midplane. The correlation between B v and 

- JB- Vu x Az is shown in the right panels of Fig. [14] with 
orange (grey) symbols indicating points above the midplane 
and blue (dark) symbols indicating points below the midplane. 
Anticorrelation below the midplane and correlation above indi- 
cates a positive a&yn above the midplane and a negative a^n be- 
low the midplane, of the order ayyn ~ 0.023c s for run S3D_256. 
The higher resolution run S3D_512 gives aa yn ~ 0.032c s . The 
fact that the helicity dynamo increases in efficiency when going 
to higher resolution, even though both the collisional hyper and 
shock resistivity coefficients are reduced, may indicate that the 
outlined dynamo is a fast dynamo, although future simulations 
applying resistivity on physical rather than on num erical grounds 
will b e needed to corroborate this point (see e.g. lHanasz et al.L 
l2002h . 

Considering the creation of radial field by the Park er insta- 
bility in a galactic environment. lHanasz & Leschl (1 1 993b predic t 
a ^ 0A{v$), where the "cyclonic velocity" (v^) of lParkerl (1 19791) 
is approximately 0.1c s in our simulations. The resulting dynamo 
coefficient o-d yn — 0.04c s is quite similar to what we find here 
from integrat ing the magnet ic stretching term. In the theoretical 
framework of Moffatt ( 1978), on the other hand, the value of ffdyn 
should be comparable to the correlation time of the turbulence 

t co1t times the mean helicity, ffd yn (1/3)t coit (V x u') ■ u', at 

least in the limit of vanishing collisional resistivity. The helicity 
in our simulations is positive below the midplane and negative 
above the midplane, with an amplitude of |(Vxm')-m'| ~ l.Thus 



the inferred sign of adyn fits well with the analytical theory, but 
the absolute value of ffdyn is at least ten times smaller than the 
expectation based on a correlation time of order unity. This dis- 
crepancy may be simply due to the fact that kinematic dynamo 
theory is not applicable to our simulations, because the Lorentz 
force plays a significant role in determining the evolution of the 
velocity field and the magnetic field. 

6. Summary and conclusions 

In this paper we consider the evolution of strongly magnetised 
Keplerian accretion discs. Our numerical experiments show that 
the hydromagnetic state of the gas flow is very different from 
what is seen in zero net flux simulations. The Parker instabil- 
ity leads to huge magnetic arcs rising several scale heights from 
the disc midplane, and the magnetorotational instability in turn 
feeds off the vertical fields and creates a highly turbulent flow, 
an inte raction that was predicted analytically by iTout & Pringld 
(119921) . Although the flow is stochastic and time fluctuating, we 
have identified an underlying dynamo process that couples the 
vertically dependent mean radial and azimuthal magnetic field 
components. As gas slides down inclined field lines, it obtains 
a helical motion due to Coriolis forces, and thus the azimuthal 
field lines are twisted in such a way as to create a mean electro- 
motoric force in the direction of the unperturbed field line - a 
configuration prone to create radial field. In turn the large scale 
radial field regenerates the azimuthal field through Keplerian 
shear. Although Parker instability dominates the linear growth 
phase, we have found evidence for magnerotational instabil- 
ity in the azimuthal field as well. In the midplane of the disc, 
where the buoyancy is weak, azimuthal MRI drives the initial 
evolution towards turbulence dFoglizzo & Tagger! fl994l fl99l 
iTerquem & Papaloizoul Il996f) . These two related instabilities, 
magnetorotational instability in the vertical fields created by the 
Parker instability and magnetorotational swing instability in the 
azimuthal fields, both rely on azimuthal flux confinement and 
can coexist in the linear as well as in the non-linear state of trans- 
magnetic accretion discs. 

Such a path to accretion, based on the interaction of Parker 
and magnetorotational instabilities, has at least two appealing 
traits. First of all that the vertical fields that feed the magne- 
torotational instability are created in a transparent way by the 
Parker instability. Zero magnetic flux models must most likely 
rely on a small scale dynamo in order to create vertical fields, 
and there is mounting evidence that such a dynamo would not 
operate in the bulk part of accretion discs where the magnetic 
Prandtl number is much lo wer than unity dSchekochihin et all 
120051 : iFromang et"aill2007l) . The second appealing result of our 
model is that the Maxwell and Reynolds stresses are significant 
(a * 0.1). Such high accretion stresses could solve the problem 
that observed accretion rates are often one or two orders of mag- 
nitude higher tha n the accretion ra tes obtained in zero net flux 
MRI simulations dKing et alll2007l) . 

The regeneration of azimuthal field by the shearing of an 
appropriate radial field was seen in all our simulations that in- 
cluded Keplerian shear. As magnetohydrostatic equilibrium is 
compromised by the Parker instability, gas streams down along 
inclined field lines. Coriolis force diverts the gas to the right, 
and a radial magnetic field is created as the azimuthal field is 
subjected to shear-regions typically the size of the Parker in- 
stability. Eventually magnetic reconnection leads to a coherent 
large s c ale rad ial magnetic field. This dynamo was predicted by 
iParkerl d 19921) and subsequently ob served in the rigid rotation 
simulations of lHanasz et al.l d2002l) . To our knowledge we are 
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the first to point out the relevance of Parker's fast galactic dy- 
namo to accretion discs and how it closes the accretion loop by 
replenishing the azimuthal field that is lost by magnetic buoy- 
ancy. 

The fine-tuned initial conditions with a purely azimuthal 
magnetic field and a constant ratio of magnetic to thermal pres- 
sure may be questioned. However our experiments with a com- 
bined azimuthal and vertical field shows that the Parker insta- 
bility is robust even if the azimuthal field coexists with a mod- 
erately strong net vertical field, and that the additional vertical 
field component may indeed increase the accretion rate further. 

Our results may also be relevant for star formation in the 
galactic centre. Although there is currently no coherent ac- 
cretion disc structure, the population of young, massive stars 
in a disc-like structure close to the galactic centre points to- 
wards the brief existence of an accretion disc some million 
years ago dLevin & Beloborodovj 2003; M ilosavlievic & Loebl 
l2004t iNavakshin et all 120071: lAlexander et all 120081) . The disc 
was likely to be initially strongly magnetised, as indicated by 
the current high magnetisation of the circumnuclear molecu- 
lar ring. Hence the type of MHD processes studied in this pa- 
per may be of central significance for the disc dynamics. The 
presented simulations of transmagnetic (J3 ~ 1) discs argues 
that discs dominated by magnetic pressure /3 <K 1 are astro- 
ph ysically viab l e. The existence of such discs was conjectured 
by iPariev et all (120031) and thei r limitations and observatio nal 
consequences were explored by iBegelman & Pringld (120071) . A 
number of problems in accretion disc theory are alleviated by 
the presence of super-equipartition magnetic fields, among them 
is the long- standing issue of self-gravity and fragmentation of 
AGN discs dGoodmanl [20031) . 

Future research into strongly magnetised accretion discs 
should also focus on the self-consistent modelling of the mag- 
netisation of the material that feeds accretion discs in various 
environments and on the evolution of magnetised disc coronae, 
in light of effects such as reconnection, shearing of foot points, 
Ohmic heating and radiative cooling that take place there. 
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